PatchSizePilot
Emanuele Giacomuzzo
2022-08-03
PatchSizePilot
Overview
Meta-ecosystems have been studied looking at meta-ecosystems in which patch size was the same. However, of course, we know that meta-ecosystems are mad out of patches that have different size. To see the effects of patch size on meta-ecosystem properties, we ran a four weeks protist experiment in which different ecosystems were connected through the flow of nutrients. The flow of nutrients resulted from a perturbation of the ecosystems in which a fixed part of the cultures was boiled and then poored into the receiving patch. This had a fixed volume (e.g., small perturbation = 6.75 ml) and was the same across all patch sizes. The experiment design consisted in crossing two disturbances with a small, medium, and large isolated ecosystems and with a small-small, medium-medium, large-large, and small-large meta-ecosystem. We took videos every four days and we create this perturbation and resource flow the day after taking videos. We skipped the perturbation the day after we assembled the experiment so that we would start perturbing it when population densities were already high.
We had mainly two research questions:
Do local properties of a patch depend upon the size of the patch it is connected to?
Do regional properties of a meta-ecosystem depend upon the relative size of its patch?
Biomass
Aim
Here I study how biomass density changes across treatments in the PatchSizePilot. In particular, my research questions are:
How does biomass density change regionally?
Do meta-ecosystems with the same total size but with patches that are either the same size or of different size have a different biomass density? (Do the medium-medium and small-large meta-ecosystems have different biomass density?)
And if they do, is it because of resource flow? Or would we see this also with small and large ecosystems that are not connected? (Do the small-large meta-ecosystems have different biomass density from two isolated small and large patches?)
- How does biomass density change locally?
- How does biomass density change according to the size to which the patch is connected? (Does a small patch connected to a small patch have the same biomass density than a small patch connected to a large patch? And does a large patch connected to a large patch have the same biomass density than a large patch connected to a small patch?)
Data
Experimental cultures
culture_info = read.csv(here("data", "PatchSizePilot_culture_info.csv"), header = TRUE)
datatable(culture_info[,1:10],
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))Local biomass dataset
### --- IMPORT --- ###
load(here("data", "population", "t0.RData")); t0 = pop_output
load(here("data", "population", "t1.RData")); t1 = pop_output
load(here("data", "population", "t2.RData")); t2 = pop_output
load(here("data", "population", "t3.RData")); t3 = pop_output
load(here("data", "population", "t4.RData")); t4 = pop_output
load(here("data", "population", "t5.RData")); t5 = pop_output
load(here("data", "population", "t6.RData")); t6 = pop_output
load(here("data", "population", "t7.RData")); t7 = pop_output
rm(pop_output)
### --- TIDY --- ###
#Column: time
t0$time = NA
t1$time = NA
#Column: replicate_video
t0$replicate_video = 1:12 #In t1 I took 12 videos of a single
t1$replicate_video = 1 #In t1 I took only 1 video/culture
t2$replicate_video = 1 #In t2 I took only 1 video/culture
t3$replicate_video = 1 #In t3 I took only 1 video/culture
t4$replicate_video = 1 #In t4 I took only 1 video/culture
t5$replicate_video = 1 #In t5 I took only 1 video/culture
t6 = t6 %>%
rename(replicate_video = replicate)
t7 = t7 %>%
rename(replicate_video = replicate)
#Create an elongated version of t0 so that each of the 110 cultures can have 12 video replicates at t0.
elongating_t0 = NULL
for (video in 1:nrow(t0)){
for (ID in 1:nrow(culture_info)) {
elongating_t0 = rbind(elongating_t0, t0[video,])
}
}
ID_vector = rep(1:nrow(culture_info),
times = nrow(t0))
elongating_t0$culture_ID = ID_vector
#Merge previous data-sets
t0 = merge(culture_info,elongating_t0, by="culture_ID")
t1 = merge(culture_info,t1, by = "culture_ID")
t2 = merge(culture_info,t2, by = "culture_ID")
t3 = merge(culture_info,t3, by = "culture_ID")
t4 = merge(culture_info,t4, by = "culture_ID")
t5 = merge(culture_info,t5, by = "culture_ID")
t6 = merge(culture_info,t6, by = "culture_ID")
t7 = merge(culture_info,t7, by = "culture_ID")
ds_biomass = rbind(t0, t1, t2, t3, t4, t5, t6, t7)
rm(elongating_t0, t0, t1, t2, t3, t4, t5, t6, t7)
#Column: time_point
ds_biomass$time_point[ds_biomass$time_point=="t0"] = 0
ds_biomass$time_point[ds_biomass$time_point=="t1"] = 1
ds_biomass$time_point[ds_biomass$time_point=="t2"] = 2
ds_biomass$time_point[ds_biomass$time_point=="t3"] = 3
ds_biomass$time_point[ds_biomass$time_point=="t4"] = 4
ds_biomass$time_point[ds_biomass$time_point=="t5"] = 5
ds_biomass$time_point[ds_biomass$time_point=="t6"] = 6
ds_biomass$time_point[ds_biomass$time_point=="t7"] = 7
ds_biomass$time_point = as.character(ds_biomass$time_point)
#System nr 40 still here
#Column: day
ds_biomass$day = NA
ds_biomass$day[ds_biomass$time_point== 0] = 0
ds_biomass$day[ds_biomass$time_point== 1] = 4
ds_biomass$day[ds_biomass$time_point== 2] = 8
ds_biomass$day[ds_biomass$time_point== 3] = 12
ds_biomass$day[ds_biomass$time_point== 4] = 16
ds_biomass$day[ds_biomass$time_point== 5] = 20
ds_biomass$day[ds_biomass$time_point== 6] = 24
ds_biomass$day[ds_biomass$time_point== 7] = 28
#Column: size_of_connected_patch
ds_biomass$size_of_connected_patch[ds_biomass$eco_metaeco_type == "S"] = "S"
ds_biomass$size_of_connected_patch[ds_biomass$eco_metaeco_type == "S (S_S)"] = "S"
ds_biomass$size_of_connected_patch[ds_biomass$eco_metaeco_type == "S (S_L)"] = "L"
ds_biomass$size_of_connected_patch[ds_biomass$eco_metaeco_type == "M (M_M)"] = "M"
ds_biomass$size_of_connected_patch[ds_biomass$eco_metaeco_type == "L"] = "L"
ds_biomass$size_of_connected_patch[ds_biomass$eco_metaeco_type == "L (L_L)"] = "L"
ds_biomass$size_of_connected_patch[ds_biomass$eco_metaeco_type == "L (S_L)"] = "S"
#Column: eco_metaeco_type
ds_biomass$eco_metaeco_type = factor(ds_biomass$eco_metaeco_type,
levels = c('S',
'S (S_S)',
'S (S_L)',
'M',
'M (M_M)',
'L',
'L (L_L)',
'L (S_L)'))
ecosystems_to_take_off = 60 #Culture number 60 because it was spilled (isolated large patch, high disturbance, system nr = 40)
ds_biomass = ds_biomass %>%
filter(! culture_ID %in% ecosystems_to_take_off)
ds_for_evaporation = ds_biomass
ds_biomass = ds_biomass %>%
select(culture_ID,
patch_size,
disturbance,
metaecosystem_type,
bioarea_per_volume,
replicate_video,
time_point,
day,
metaecosystem,
system_nr,
eco_metaeco_type,
size_of_connected_patch) %>%
relocate(culture_ID,
system_nr,
disturbance,
time_point,
day,
patch_size,
metaecosystem,
metaecosystem_type,
eco_metaeco_type,
size_of_connected_patch,
replicate_video,
bioarea_per_volume)
datatable(ds_biomass,
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))Regional biomass data set
ds_regional = ds_biomass %>%
filter(metaecosystem == "yes") %>%
group_by(culture_ID,
system_nr,
disturbance,
day,
time_point,
patch_size,
metaecosystem_type) %>%
summarise(patch_mean_bioarea_across_videos = mean(bioarea_per_volume)) %>%
group_by(system_nr, disturbance, day, time_point, metaecosystem_type) %>%
summarise(regional_mean_bioarea = mean(patch_mean_bioarea_across_videos))
metaecosystems_to_take_off = 40 #System 40 was the system of culture 60 that I spilled
ds_regional = ds_regional %>%
filter(! system_nr %in% metaecosystems_to_take_off)
datatable(ds_regional,
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))Local biomass lnRR data-set
for (disturbance_input in c("low", "high")){
for (eco_metaeco_input in c("S", "M", "L")){
for (time_point_input in 0:7){
ds_biomass$isolated_control[ds_biomass$patch_size == eco_metaeco_input &
ds_biomass$disturbance == disturbance_input &
ds_biomass$time_point == time_point_input] =
ds_biomass %>%
filter(disturbance == disturbance_input) %>%
filter(eco_metaeco_type == eco_metaeco_input) %>%
filter(time_point == time_point_input) %>%
group_by(culture_ID) %>%
summarise(bioarea_per_volume_across_videos = mean(bioarea_per_volume)) %>%
summarise(mean_bioarea_per_volume = mean(bioarea_per_volume_across_videos))
}
}
}
ds_biomass_averaged_across_videos = ds_biomass %>%
mutate(isolated_control = as.numeric(isolated_control)) %>%
group_by(disturbance,
eco_metaeco_type,
culture_ID,
day,
time_point,
isolated_control) %>%
summarise(bioarea_per_volume_across_videos = mean(bioarea_per_volume)) %>% #Across videos
mutate(RR_bioarea_per_volume = (bioarea_per_volume_across_videos+1) /isolated_control) %>%
mutate(lnRR_bioarea_per_volume = ln(RR_bioarea_per_volume))
datatable(ds_biomass_averaged_across_videos,
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))Regional biomass
Medium-Medium vs Small-Large
Do meta-ecosystems with the same total size but with patches that are either the same size or of different size have a different biomass density? (Do the medium-medium and small-large meta-ecosystems have different biomass density?)
Plots
ds_regional %>%
filter ( disturbance == "low") %>%
filter (metaecosystem_type == "S_L" |
metaecosystem_type == "M_M") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
fill = "System nr",
color = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_linetype_discrete(labels = c("medium-medium",
"small-large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional %>%
filter ( disturbance == "high") %>%
filter (metaecosystem_type == "S_L" | metaecosystem_type == "M_M") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = high",
fill = "System nr",
color = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_linetype_discrete(labels = c("medium-medium",
"small-large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
regional_publication = ds_regional %>%
filter(disturbance == "low") %>%
filter (metaecosystem_type == "S_L" | metaecosystem_type == "M_M") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
color='',
fill='') +
scale_fill_discrete(labels = c("medium-medium",
"small-large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
regional_publication
ds_regional %>%
filter(disturbance == "high") %>%
filter (metaecosystem_type == "S_L" | metaecosystem_type == "M_M") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = interaction (day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = high",
color='',
fill='') +
scale_fill_discrete(labels = c("medium-medium", "small-large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
Model time series
How does the biomass density of medium-medium and small-large meta-ecosystems differ across the time series? (The first two points before the first disturbance are taken off).
Let’s see how linear is the time trend of bioarea and if we can make it more linear with a log10 transformation. We are lucky that during the modelling process we need to drop the first two time points because they were before the first perturbation.
Linearity of regional bioarea ~ time
ds_regional %>%
filter(time_point >= 2) %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = day)) +
geom_boxplot() +
labs(title = "Without log transformation",
x = "Day",
y = "Regional bioarea (something/µl)")
linear_model = lm(regional_mean_bioarea ~
day,
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"))
par(mfrow=c(2,3))
plot(linear_model, which = 1:5)
When we don’t transform it, the time trend of bioarea is not that linear, as we can see in the problematic banana shape of the line in the residuals vs fitted plot.
Linearity of Log10(Regional bioarea +1) ~ time
ds_regional %>%
filter(time_point >= 2) %>%
ggplot(aes(x = day,
y = log(regional_mean_bioarea + 1),
group = day)) +
geom_boxplot() +
labs(title = "With log transformation",
x = "Day",
y = "Log (regional bioarea + 1) (something/µl)")
log_linear_model = lm(log10(regional_mean_bioarea + 1) ~
day,
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"))
par(mfrow=c(2,3))
plot(log_linear_model, which = 1:5)
par(mfrow=c(1,1))
The log transformed bioarea is linear.
Model selection
Let’ start from the full model.
\[ Log_{10} (Regional \: bioarea + 1) = t + M + D + tM + tD + MD + tDM + (t | system \: nr) \]
full = lmer(log10(regional_mean_bioarea + 1) ~
day * metaecosystem_type * disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))Should we keep the correlation in (day | system_nr)?
no_correlation = lmer(log10(regional_mean_bioarea + 1) ~
day * metaecosystem_type * disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(full, no_correlation)## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## full: log10(regional_mean_bioarea + 1) ~ day * metaecosystem_type * disturbance + (day | system_nr)
## no_correlation: log10(regional_mean_bioarea + 1) ~ day * metaecosystem_type * disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## full 12 -168.78 -135.33 96.389 -192.78
## no_correlation 12 -168.78 -135.33 96.389 -192.78 0 0Yes.
Should we keep t * M * D?
no_threeway = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : metaecosystem_type +
day : disturbance +
metaecosystem_type : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = 'optimx',
optCtrl = list(method = 'L-BFGS-B')))
anova(full, no_threeway)## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_threeway: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:metaecosystem_type + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## full: log10(regional_mean_bioarea + 1) ~ day * metaecosystem_type * disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_threeway 11 -170.66 -140.00 96.330 -192.66
## full 12 -168.78 -135.33 96.389 -192.78 0.1182 1 0.7309No.
Should we keep t * M?
no_TM = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : disturbance +
metaecosystem_type : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(no_threeway,no_TM)## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_TM: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## no_threeway: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:metaecosystem_type + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_TM 10 -172.62 -144.74 96.308 -192.62
## no_threeway 11 -170.66 -140.00 96.330 -192.66 0.0431 1 0.8356No.
Should we keep t * D?
no_TD = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
metaecosystem_type : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(no_TM, no_TD)## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_TD: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + metaecosystem_type:disturbance + (day | system_nr)
## no_TM: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_TD 9 -163.54 -138.46 90.771 -181.54
## no_TM 10 -172.62 -144.74 96.308 -192.62 11.074 1 0.0008756 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1No.
Should we keep M * D?
no_MD = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(no_TM, no_MD)## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_MD: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + (day | system_nr)
## no_TM: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + metaecosystem_type:disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_MD 9 -173.87 -148.78 95.932 -191.87
## no_TM 10 -172.62 -144.74 96.308 -192.62 0.7513 1 0.3861No.
Should we keep the random effect of system nr on the time slopes (day | system_nr)?
no_random_slopes = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : disturbance +
(1 | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
anova(no_MD, no_random_slopes)## Data: ds_regional %>% filter(time_point >= 2) %>% filter(metaecosystem_type == ...
## Models:
## no_random_slopes: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + (1 | system_nr)
## no_MD: log10(regional_mean_bioarea + 1) ~ day + metaecosystem_type + disturbance + day:disturbance + (day | system_nr)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_random_slopes 7 -173.23 -153.71 93.613 -187.23
## no_MD 9 -173.87 -148.78 95.932 -191.87 4.6394 2 0.0983 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Yes.
Best model
Therefore, our best model is:
\[ log_{10}(regional \: bioarea + 1) = t + M + D + t*D + (t| system \: nr) \]
The R squared of this model for t2-t7 are:
best_model = no_MD
R2_marginal = r.squaredGLMM(best_model)[1]
R2_marginal = round(R2_marginal, digits = 2)
R2_conditional = r.squaredGLMM(best_model)[2]
R2_conditional = round(R2_conditional, digits = 2)Marginal R2 = 0.84
Conditional R2 = 0.88
Let’s just assume that this model holds also for t2-t5. Then, let’s recalculate the R squared.
t2_t5 = lmer(log10(regional_mean_bioarea + 1) ~
day +
metaecosystem_type +
disturbance +
day : disturbance +
(day | system_nr),
data = ds_regional %>%
filter(time_point >= 2) %>%
filter(time_point <= 5) %>%
filter(metaecosystem_type == "M_M" | metaecosystem_type == "S_L"),
REML = FALSE,
control = lmerControl(optimizer = "Nelder_Mead"))
R2_marginal = r.squaredGLMM(t2_t5)[1]
R2_marginal = round(R2_marginal, digits = 2)
R2_conditional = r.squaredGLMM(t2_t5)[2]
R2_conditional = round(R2_conditional, digits = 2)The R squared of this model for t2-t5 are:
Marginal R2 = 0.7
Conditional R2 = 0.8
Next steps:
- calculating the R2 of the meta-ecosystem type. I was thinking I
could do it with
part2but it can’t do it when the model has a random slope. -> evaluating whether meta-ecosystem type R2 is high enough -> if it’s not, redo model selection with t2-t5 to see if you can increase R2 of meta-ecosystem type.
Model single points
How does the biomass density of medium-medium and small-large meta-ecosystems differ for each time point? (The first two points before the first disturbance are taken off).
Need to wait until I understand how to calculate R2 for models with random slopes, so then I can look at the R2 of meta-ecosystem type across time points.
full = lmer(log10(regional_mean_bioarea + 1) ~
metaecosystem_type +
disturbance +
metaecosystem_type : disturbance +
(1 | system_nr),
data = ds_regional %>%
filter(time_point == 3) %>%
filter(metaecosystem_type == "M_M" |
metaecosystem_type == "S_L"),
REML = FALSE)
# Next step: understanding what is a grouping factor. I need to understand what a grouping factor because I cannot run the full model, as it gives me the following problem: "Error: number of levels of each grouping factor must be < number of observations (problems: system_nr)."Small-Large vs Small-Large from isolated
Do a meta-ecosystem with patches of the same size and a meta-ecosystems with patches of different size have different regional biomass density because of their resource flow? Or would we see this also with small and large ecosystems that are not connected? (Do the small-large meta-ecosystems have different biomass density from two isolated small and large patches?)
Tidy
To create the combinations of small and large isolated patches to be compared to the small-large meta-ecosystems, I’m going through the following steps.
- Take isolated small and large patches.
- Because I had to take out the last
- Add them together randomly (here following numeric order, e.g., small isolated 1 with large isolated 11, small isolated 2 with large isolated 12 and so on). Doing this, we make sure that we are matching their disturbances.
system_nr_to_take_off = 50
n_time_points = 8
isolated_S_and_L = ds_biomass %>%
filter(!system_nr == system_nr_to_take_off) %>%
filter(eco_metaeco_type == "S" | eco_metaeco_type == "L") %>%
group_by(system_nr, disturbance, time_point, day, eco_metaeco_type) %>%
summarise(bioarea_per_volume_across_videos = mean(bioarea_per_volume))
### --- Low disturbance --- ###
isolated_S_low = isolated_S_and_L %>%
filter(eco_metaeco_type == "S") %>%
filter(disturbance == "low")
isolated_L_low = isolated_S_and_L %>%
filter(eco_metaeco_type == "L") %>%
filter(disturbance == "low")
n_isolated_patches = (nrow(isolated_S_low)) / n_time_points
number_for_pairing_v = rep(1:n_isolated_patches, each = n_time_points)
isolated_S_low$number_for_pairing = number_for_pairing_v
isolated_L_low$number_for_pairing = number_for_pairing_v
### --- High disturbance --- ###
isolated_S_high = isolated_S_and_L %>%
filter(eco_metaeco_type == "S") %>%
filter(disturbance == "high")
isolated_L_high = isolated_S_and_L %>%
filter(eco_metaeco_type == "L") %>%
filter(disturbance == "high")
n_isolated_patches = (nrow(isolated_S_high)) / n_time_points
number_for_pairing_v = rep(1:n_isolated_patches, each = n_time_points)
isolated_S_high$number_for_pairing = number_for_pairing_v + 100
isolated_L_high$number_for_pairing = number_for_pairing_v + 100
SL_from_isolated = rbind(isolated_S_low,
isolated_L_low,
isolated_S_high,
isolated_L_high) %>%
mutate(metaecosystem_type = "S_L_from_isolated") %>%
rename(regional_mean_bioarea = bioarea_per_volume_across_videos)
ds_regional_with_SL_from_isolated = rbind(SL_from_isolated, ds_regional)#This would be another way of doing it
system_nr_S_low = unique(isolated_S$system_nr)[1:5]
system_nr_L_low = unique(isolated_L$system_nr)[1:5]
system_nr_S_high = unique(isolated_S$system_nr)[6:9]
system_nr_L_high = unique(isolated_L$system_nr)[6:9]
low_pairs = expand.grid(system_nr_S_low,system_nr_L_low)
high_pairs = expand.grid(system_nr_S_high, system_nr_L_high)
pairs = rbind(low_pairs, high_pairs)
number_of_pairs = nrow(pairs)
SL_from_isolated_all_combinations = NULL
for (pair in 1:number_of_pairs){
SL_from_isolated_one_combination = ds_biomass %>%
filter(system_nr %in% pairs[pair,]) %>%
group_by(disturbance, day, time_point, system_nr) %>%
summarise(regional_bioarea_across_videos = mean(bioarea_per_volume)) %>%
group_by(disturbance, day, time_point) %>%
summarise(regional_mean_bioarea = mean(regional_bioarea_across_videos)) %>%
mutate(system_nr = 1000 + pair) %>%
mutate(metaecosystem_type = "S_L_from_isolated")
SL_from_isolated_all_combinations = rbind(SL_from_isolated_one_combination,
SL_from_isolated_all_combinations)
}
ds_regional_with_SL_from_isolated = rbind(SL_from_isolated_all_combinations, ds_regional)Next step: bootstrap 5 isolated small patches and 5 isolated large patches
Plots
ds_regional_with_SL_from_isolated %>%
filter ( disturbance == "low") %>%
filter (metaecosystem_type == "S_L" | metaecosystem_type == "S_L_from_isolated") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
fill = "System nr",
color = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_linetype_discrete(labels = c("small-large",
"small-large \n from isolated")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional_with_SL_from_isolated %>%
filter ( disturbance == "high") %>%
filter (metaecosystem_type == "S_L" | metaecosystem_type == "S_L_from_isolated") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(title = "Disturbance = high",
x = "Day",
y = "Regional bioarea (something/µl)",
fill = "System nr",
color = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_linetype_discrete(labels = c("small-large",
"small-large \n from isolated")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional_with_SL_from_isolated %>%
filter(disturbance == "low") %>%
filter(metaecosystem_type == "S_L" | metaecosystem_type == "S_L_from_isolated") %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(title = "Disturbance = low",
x = "Day",
y = "Regional bioarea (something/microlitre)",
fill = "") +
scale_fill_discrete(labels = c("small-large", "isolated small & \n isolated large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional_with_SL_from_isolated %>%
filter(disturbance == "high") %>%
filter(metaecosystem_type == "S_L" | metaecosystem_type == "S_L_from_isolated") %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "Regional bioarea (something/microlitre)",
fill = "") +
scale_fill_discrete(labels = c("small-large", "isolated small & \n isolated large")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")The plots broke when I tried to do the bootstrapping.
Modelling
The model broke when I tried to do the bootstrapping.
mixed_model = lmer(log10(regional_mean_bioarea +1) ~
day * metaecosystem_type * disturbance +
(day | system_nr),
data = ds_regional_with_SL_from_isolated %>%
filter(metaecosystem_type == "S_L" |
metaecosystem_type == "S_L_from_isolated") %>%
filter(time_point >= 2),
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
null_model = lmer(log10(regional_mean_bioarea +1) ~
day * disturbance +
(day | system_nr),
data = ds_regional_with_SL_from_isolated %>%
filter(metaecosystem_type == "S_L" |
metaecosystem_type == "S_L_from_isolated") %>%
filter(time_point >= 2),
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
anova(mixed_model, null_model)Meta-ecosystems of different total size
How does the biomass density of meta-ecosystems change according to the size of their patches?
ds_regional %>%
filter(!metaecosystem_type == "S_L") %>%
filter ( disturbance == "low") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
fill = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_colour_continuous(guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large-large",
"medium-medium",
"small-small")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional %>%
filter(!metaecosystem_type == "S_L") %>%
filter ( disturbance == "high") %>%
ggplot (aes(x = day,
y = regional_mean_bioarea,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = metaecosystem_type)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = high",
fill = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_colour_continuous(guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large-large",
"medium-medium",
"small-small")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional %>%
filter(disturbance == "low") %>%
filter(!metaecosystem_type == "S_L") %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(title = "Disturbance = low",
x = "Day",
y = "Regional bioarea (something/μl)",
fill = "") +
#scale_fill_discrete(labels = c("isolated large", "isolated medium", "isolated small")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large-large",
"medium-medium",
"small-small")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_regional %>%
filter(disturbance == "high") %>%
filter(!metaecosystem_type == "S_L") %>%
ggplot(aes(x = day,
y = regional_mean_bioarea,
group = interaction(day, metaecosystem_type),
fill = metaecosystem_type)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "Regional bioarea (something/μl)",
fill = "") +
#scale_fill_discrete(labels = c("isolated large", "isolated medium", "isolated small")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large-large",
"medium-medium",
"small-small")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
Interesting. It seems like there’s not much difference between the medium-medium and the large-large.
Local biomass
Small patches
How does biomass density change according to the size to which the patch is connected? (Does a small patch connected to a small patch have the same biomass density than a small patch connected to a large patch?)
Plots
ds_biomass %>%
filter(disturbance == "low") %>%
filter(patch_size == "S") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = culture_ID,
fill = culture_ID,
color = culture_ID,
linetype = eco_metaeco_type)) +
geom_line(stat = "summary", fun = "mean") +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = low",
linetype = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("small isolated",
"small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_biomass %>%
filter(disturbance == "high") %>%
filter(patch_size == "S") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = culture_ID,
fill = system_nr,
color = system_nr,
linetype = eco_metaeco_type)) +
geom_line(stat = "summary", fun = "mean") +
labs(title = "Disturbance = high",
x = "Day",
y = "Local bioarea (something/μl)",
linetype = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("small isolated",
"small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_biomass %>%
filter(disturbance == "low") %>%
filter(patch_size == "S") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day,eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = low",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("small isolated",
"small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_biomass %>%
filter(disturbance == "high") %>%
filter(patch_size == "S") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day,eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "Local bioarea (something/μl)",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("small isolated",
"small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
small_patches_publication = ds_biomass_averaged_across_videos %>%
filter(disturbance == "low") %>%
filter(eco_metaeco_type == "S (S_L)" | eco_metaeco_type == "S (S_S)") %>%
ggplot(aes(x = day,
y = lnRR_bioarea_per_volume,
group = interaction(day, eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(title = "Disturbance = low",
x = "Day",
y = "LnRR Local Bioarea (something/μl)",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.3, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation, LnRR Local Bioarea: Bioarea + 1 / Mean bioarea small isolated")
small_patches_publication
ds_biomass_averaged_across_videos %>%
filter(disturbance == "high") %>%
filter(eco_metaeco_type == "S (S_L)" | eco_metaeco_type == "S (S_S)") %>%
ggplot(aes(x = day,
y = lnRR_bioarea_per_volume,
group = interaction(day, eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "LnRR Local Bioarea (something/μl)",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.3, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("small connected to small",
"small connected to large")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation, LnRR Local Bioarea: Bioarea/Mean bioarea small isolated")
Time series
Model selection
Let’s start from the full model:
\[ ln \: RR (bioarea) = t + M + D + t*M + t*D + M*D + t*M+D + (t | system \: nr) \]
And let’s consider all data points from the second to the seventh.
first_time_point = 2
last_time_point = 7full_model = lmer(lnRR_bioarea_per_volume ~
day * eco_metaeco_type * disturbance +
(day | culture_ID),
data = ds_biomass_averaged_across_videos %>%
filter(time_point >= first_time_point) %>%
filter(time_point <= last_time_point) %>%
filter(eco_metaeco_type== "S (S_S)" |
eco_metaeco_type == "S (S_L)"),
REML = FALSE)Should we keep the interaction between the intercept and slope in (day | culture_ID)?
no_interaction = lmer(lnRR_bioarea_per_volume ~
day * eco_metaeco_type * disturbance +
(day || culture_ID),
data = ds_biomass_averaged_across_videos %>%
filter(time_point >= first_time_point) %>%
filter(time_point <= last_time_point) %>%
filter(eco_metaeco_type== "S (S_S)" |
eco_metaeco_type == "S (S_L)"),
REML = FALSE)
anova(full_model, no_interaction)## Data: ds_biomass_averaged_across_videos %>% filter(time_point >= first_time_point) %>% ...
## Models:
## no_interaction: lnRR_bioarea_per_volume ~ day * eco_metaeco_type * disturbance + ((1 | culture_ID) + (0 + day | culture_ID))
## full_model: lnRR_bioarea_per_volume ~ day * eco_metaeco_type * disturbance + (day | culture_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_interaction 11 465.78 500.90 -221.89 443.78
## full_model 12 461.92 500.23 -218.96 437.92 5.8649 1 0.01545 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Yes.
Should we keep the effects on the slope?
no_slope = lmer(lnRR_bioarea_per_volume ~
day * eco_metaeco_type * disturbance +
(1 | culture_ID),
data = ds_biomass_averaged_across_videos %>%
filter(time_point >= first_time_point) %>%
filter(time_point <= last_time_point) %>%
filter(eco_metaeco_type== "S (S_S)" |
eco_metaeco_type == "S (S_L)"),
REML = FALSE)
anova(full_model, no_slope)## Data: ds_biomass_averaged_across_videos %>% filter(time_point >= first_time_point) %>% ...
## Models:
## no_slope: lnRR_bioarea_per_volume ~ day * eco_metaeco_type * disturbance + (1 | culture_ID)
## full_model: lnRR_bioarea_per_volume ~ day * eco_metaeco_type * disturbance + (day | culture_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_slope 10 472.17 504.10 -226.08 452.17
## full_model 12 461.92 500.23 -218.96 437.92 14.248 2 0.0008055 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Yes.
Should we keep t * M * D?
no_three_way = lmer(lnRR_bioarea_per_volume ~
day +
eco_metaeco_type +
disturbance +
day : eco_metaeco_type +
day : disturbance +
eco_metaeco_type : disturbance +
(day | culture_ID),
data = ds_biomass_averaged_across_videos %>%
filter(time_point >= first_time_point) %>%
filter(time_point <= last_time_point) %>%
filter(eco_metaeco_type== "S (S_S)" |
eco_metaeco_type == "S (S_L)"),
REML = FALSE)
anova(full_model, no_three_way)## Data: ds_biomass_averaged_across_videos %>% filter(time_point >= first_time_point) %>% ...
## Models:
## no_three_way: lnRR_bioarea_per_volume ~ day + eco_metaeco_type + disturbance + day:eco_metaeco_type + day:disturbance + eco_metaeco_type:disturbance + (day | culture_ID)
## full_model: lnRR_bioarea_per_volume ~ day * eco_metaeco_type * disturbance + (day | culture_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_three_way 11 460.13 495.25 -219.06 438.13
## full_model 12 461.92 500.23 -218.96 437.92 0.2107 1 0.6462No.
Should we keep t * M?
no_TM = lmer(lnRR_bioarea_per_volume ~
day +
eco_metaeco_type +
disturbance +
day : disturbance +
eco_metaeco_type : disturbance +
(day | culture_ID),
data = ds_biomass_averaged_across_videos %>%
filter(time_point >= first_time_point) %>%
filter(time_point <= last_time_point) %>%
filter(eco_metaeco_type== "S (S_S)" |
eco_metaeco_type == "S (S_L)"),
REML = FALSE)
anova(no_three_way, no_TM)## Data: ds_biomass_averaged_across_videos %>% filter(time_point >= first_time_point) %>% ...
## Models:
## no_TM: lnRR_bioarea_per_volume ~ day + eco_metaeco_type + disturbance + day:disturbance + eco_metaeco_type:disturbance + (day | culture_ID)
## no_three_way: lnRR_bioarea_per_volume ~ day + eco_metaeco_type + disturbance + day:eco_metaeco_type + day:disturbance + eco_metaeco_type:disturbance + (day | culture_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_TM 10 466.49 498.42 -223.25 446.49
## no_three_way 11 460.13 495.25 -219.06 438.13 8.3612 1 0.003833 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Yes.
Should we keep t * D?
no_TD = lmer(lnRR_bioarea_per_volume ~
day +
eco_metaeco_type +
disturbance +
day : eco_metaeco_type +
eco_metaeco_type : disturbance +
(day | culture_ID),
data = ds_biomass_averaged_across_videos %>%
filter(time_point >= first_time_point) %>%
filter(time_point <= last_time_point) %>%
filter(eco_metaeco_type== "S (S_S)" |
eco_metaeco_type == "S (S_L)"),
REML = FALSE)
anova(no_three_way, no_TD)## Data: ds_biomass_averaged_across_videos %>% filter(time_point >= first_time_point) %>% ...
## Models:
## no_TD: lnRR_bioarea_per_volume ~ day + eco_metaeco_type + disturbance + day:eco_metaeco_type + eco_metaeco_type:disturbance + (day | culture_ID)
## no_three_way: lnRR_bioarea_per_volume ~ day + eco_metaeco_type + disturbance + day:eco_metaeco_type + day:disturbance + eco_metaeco_type:disturbance + (day | culture_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_TD 10 458.86 490.79 -219.43 438.86
## no_three_way 11 460.13 495.25 -219.06 438.13 0.7309 1 0.3926No.
Should we keep M * D?
no_MD = lmer(lnRR_bioarea_per_volume ~
day +
eco_metaeco_type +
disturbance +
day : eco_metaeco_type +
(day | culture_ID),
data = ds_biomass_averaged_across_videos %>%
filter(time_point >= first_time_point) %>%
filter(time_point <= last_time_point) %>%
filter(eco_metaeco_type== "S (S_S)" |
eco_metaeco_type == "S (S_L)"),
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
anova(no_TD, no_MD)## Data: ds_biomass_averaged_across_videos %>% filter(time_point >= first_time_point) %>% ...
## Models:
## no_MD: lnRR_bioarea_per_volume ~ day + eco_metaeco_type + disturbance + day:eco_metaeco_type + (day | culture_ID)
## no_TD: lnRR_bioarea_per_volume ~ day + eco_metaeco_type + disturbance + day:eco_metaeco_type + eco_metaeco_type:disturbance + (day | culture_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_MD 9 456.94 485.68 -219.47 438.94
## no_TD 10 458.86 490.79 -219.43 438.86 0.0842 1 0.7717No.
Should we keep D?
no_D = lmer(lnRR_bioarea_per_volume ~
day +
eco_metaeco_type +
day : eco_metaeco_type +
(day | culture_ID),
data = ds_biomass_averaged_across_videos %>%
filter(time_point >= first_time_point) %>%
filter(time_point <= last_time_point) %>%
filter(eco_metaeco_type== "S (S_S)" |
eco_metaeco_type == "S (S_L)"),
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead")
)
anova(no_MD, no_D)## Data: ds_biomass_averaged_across_videos %>% filter(time_point >= first_time_point) %>% ...
## Models:
## no_D: lnRR_bioarea_per_volume ~ day + eco_metaeco_type + day:eco_metaeco_type + (day | culture_ID)
## no_MD: lnRR_bioarea_per_volume ~ day + eco_metaeco_type + disturbance + day:eco_metaeco_type + (day | culture_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_D 8 466.08 491.62 -225.04 450.08
## no_MD 9 456.94 485.68 -219.47 438.94 11.137 1 0.0008464 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Yes.
Best model
Then our best model is:
\[ RR (bioarea) = t + M + D + t*M + (t | ID) \]
Let’s then do some model diagnostics.
best_model = no_MD
plot(best_model)
qqnorm(resid(best_model))
The R squared of this model for t2-t7 are:
R2_marginal = r.squaredGLMM(best_model)[1]
R2_marginal = round(R2_marginal, digits = 2)
R2_conditional = r.squaredGLMM(best_model)[2]
R2_conditional = round(R2_conditional, digits = 2)Marginal R2 = 0.22
Conditional R2 = 0.4
Let’s just assume that this model holds also for t2-t5. Then, let’s recalculate the R squared.
last_time_point = 5
t2_t5 = lmer(lnRR_bioarea_per_volume ~
day +
eco_metaeco_type +
disturbance +
day : eco_metaeco_type +
(day | culture_ID),
data = ds_biomass_averaged_across_videos %>%
filter(time_point >= first_time_point) %>%
filter(time_point <= last_time_point) %>%
filter(eco_metaeco_type== "S (S_S)" |
eco_metaeco_type == "S (S_L)"),
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))The R squared of this model for t2-t5 are:
R2_marginal = r.squaredGLMM(t2_t5)[1]
R2_marginal = round(R2_marginal, digits = 2)
R2_conditional = r.squaredGLMM(t2_t5)[2]
R2_conditional = round(R2_conditional, digits = 2)Marginal R2 = 0.38
Conditional R2 = 0.46
Single time points
- Next step: understanding how to calculate the R2 of meta-ecosystem type when there are random slopes.
Large patches
How does biomass density change according to the size to which the patch is connected? (Does a large patch connected to a large patch have the same biomass density than a large patch connected to a small patch?)
plots
ds_biomass %>%
filter(disturbance == "low") %>%
filter(patch_size == "L") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = eco_metaeco_type)) +
geom_line(stat = "summary", fun = "mean") +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = low",
linetype = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large isolated",
"large connected to large",
"large connected to small")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_biomass %>%
filter(disturbance == "high") %>%
filter(patch_size == "L") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = eco_metaeco_type)) +
geom_line(stat = "summary", fun = "mean") +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = high",
linetype = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large isolated",
"large connected to large",
"large connected to small")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
ds_biomass %>%
filter(disturbance == "low") %>%
filter(patch_size == "L") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day,eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = low",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large isolated",
"large connected to large",
"large connected to small")) +
geom_vline(xintercept = first_perturbation_day,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
local_large_high_plot = ds_biomass %>%
filter(disturbance == "high") %>%
filter(patch_size == "L") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day,eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(x = "Day",
y = "Local bioarea (something/μl)",
title = "Disturbance = high",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large isolated",
"large connected to large",
"large connected to small")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation")
local_large_high_plot
large_patches_publication = ds_biomass_averaged_across_videos %>%
filter(disturbance == "low") %>%
filter(eco_metaeco_type == "L (L_L)" | eco_metaeco_type == "L (S_L)") %>%
ggplot(aes(x = day,
y = lnRR_bioarea_per_volume,
group = interaction(day, eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(title = "Disturbance = low",
x = "Day",
y = "LnRR Local Bioarea (something/μl)",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.9, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large connected to large",
"large connected to small")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation,
Response ratio bioarea: Bioarea/Mean bioarea small isolated")
large_patches_publication
ds_biomass_averaged_across_videos %>%
filter(disturbance == "high") %>%
filter(eco_metaeco_type == "L (L_L)" | eco_metaeco_type == "L (S_L)") %>%
ggplot(aes(x = day,
y = lnRR_bioarea_per_volume,
group = interaction(day, eco_metaeco_type),
fill = eco_metaeco_type)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "LnRR Local Bioarea (something/μl)",
fill = "") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.9, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_fill_discrete(labels = c("large connected to large",
"large connected to small")) +
geom_vline(xintercept = first_perturbation_day + 0.7,
linetype="dotdash",
color = "grey",
size=0.7) +
labs(caption = "Vertical grey line: first perturbation,
Response ratio bioarea: Bioarea/Mean bioarea small isolated")
isolated vs connected to large
Does M have an effect?
full_model = lmer(bioarea_per_volume ~
metaecosystem_type * disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L" |
eco_metaeco_type == "L (L_L)") %>%
filter(time_point >= 2),
REML = FALSE)
no_metaeco_type_model = lmer(bioarea_per_volume ~
disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L" |
eco_metaeco_type == "L (L_L)") %>%
filter(time_point >= 2),
REML = FALSE)
anova(full_model, no_metaeco_type_model)## Data: ds_biomass %>% filter(eco_metaeco_type == "L" | eco_metaeco_type == ...
## Models:
## no_metaeco_type_model: bioarea_per_volume ~ disturbance + (1 | system_nr) + (1 | day)
## full_model: bioarea_per_volume ~ metaecosystem_type * disturbance + (1 | system_nr) + (1 | day)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_metaeco_type_model 5 3879.5 3896.9 -1934.7 3869.5
## full_model 7 3876.5 3900.9 -1931.2 3862.5 6.9701 2 0.03065
##
## no_metaeco_type_model
## full_model *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Yes.
isolated vs connected to small
Does M have an effect?
full_model = lmer(bioarea_per_volume ~
metaecosystem_type * disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L" |
eco_metaeco_type == "L (S_L)") %>%
filter(time_point >= 2),
REML = FALSE)
no_metaeco_type_model = lmer(bioarea_per_volume ~
disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L" |
eco_metaeco_type == "L (S_L)") %>%
filter(time_point >= 2),
REML = FALSE)
anova(full_model, no_metaeco_type_model)## Data: ds_biomass %>% filter(eco_metaeco_type == "L" | eco_metaeco_type == ...
## Models:
## no_metaeco_type_model: bioarea_per_volume ~ disturbance + (1 | system_nr) + (1 | day)
## full_model: bioarea_per_volume ~ metaecosystem_type * disturbance + (1 | system_nr) + (1 | day)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_metaeco_type_model 5 2566.9 2582.2 -1278.4 2556.9
## full_model 7 2559.7 2581.2 -1272.9 2545.7 11.151 2 0.00379
##
## no_metaeco_type_model
## full_model **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Yes.
connected to large vs connected to small
Does M have an effect?
full_model = lmer(bioarea_per_volume ~
metaecosystem_type * disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L (L_L)" |
eco_metaeco_type == "L (S_L)") %>%
filter(time_point >= 2),
REML = FALSE)
no_metaeco_type_model = lmer(bioarea_per_volume ~
disturbance +
(1 | system_nr) +
(1 | day),
data = ds_biomass %>%
filter (eco_metaeco_type == "L (L_L)" |
eco_metaeco_type == "L (S_L)") %>%
filter(time_point >= 2),
REML = FALSE)
anova(full_model, no_metaeco_type_model)## Data: ds_biomass %>% filter(eco_metaeco_type == "L (L_L)" | eco_metaeco_type == ...
## Models:
## no_metaeco_type_model: bioarea_per_volume ~ disturbance + (1 | system_nr) + (1 | day)
## full_model: bioarea_per_volume ~ metaecosystem_type * disturbance + (1 | system_nr) + (1 | day)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## no_metaeco_type_model 5 3835.5 3852.9 -1912.8 3825.5
## full_model 7 3832.2 3856.6 -1909.1 3818.2 7.2508 2 0.02664
##
## no_metaeco_type_model
## full_model *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Yes.
Isolated patches
How does biomass density change according to the size of isolated patches? (How does the biomass of small, medium, and large patches change?)
ds_biomass %>%
filter ( disturbance == "low") %>%
filter(metaecosystem == "no") %>%
group_by (system_nr, day, patch_size) %>%
summarise(mean_bioarea_per_volume_across_videos = mean(bioarea_per_volume)) %>%
ggplot (aes(x = day,
y = mean_bioarea_per_volume_across_videos,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = patch_size)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
fill = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_colour_continuous(guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large isolated",
"medium isolated",
"small isolated"))
ds_biomass %>%
filter ( disturbance == "high") %>%
filter(metaecosystem == "no") %>%
group_by (system_nr, day, patch_size) %>%
summarise(mean_bioarea_per_volume_across_videos = mean(bioarea_per_volume)) %>%
ggplot (aes(x = day,
y = mean_bioarea_per_volume_across_videos,
group = system_nr,
fill = system_nr,
color = system_nr,
linetype = patch_size)) +
geom_line () +
labs(x = "Day",
y = "Regional bioarea (something/µl)",
title = "Disturbance = low",
fill = "System nr",
linetype = "") +
scale_y_continuous(limits = c(0, 6250)) +
scale_x_continuous(limits = c(-2, 30)) +
scale_colour_continuous(guide = "none") +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
scale_linetype_discrete(labels = c("large isolated",
"medium isolated",
"small isolated"))
ds_biomass %>%
filter(disturbance == "low") %>%
filter(metaecosystem == "no") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day, patch_size),
fill = patch_size)) +
geom_boxplot() +
labs(title = "Disturbance = low",
x = "Day",
y = "Local bioarea (something/μl)",
fill = "") +
scale_fill_discrete(labels = c("isolated large", "isolated medium", "isolated small")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6))
ds_biomass %>%
filter(disturbance == "high") %>%
filter(metaecosystem == "no") %>%
ggplot(aes(x = day,
y = bioarea_per_volume,
group = interaction(day, patch_size),
fill = patch_size)) +
geom_boxplot() +
labs(title = "Disturbance = high",
x = "Day",
y = "Local bioarea (something/μl)",
fill = "") +
scale_fill_discrete(labels = c("isolated large", "isolated medium", "isolated small")) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6))
Evaporation
Was the amount of evaporation different across different treatments and time points?
We want to know if there was a systematic bias in the evaporation of different treatments (disturbance, patch size) and whether evaporation changed across time. My expectation would be that we would see a difference among the exchanges 2,3 and the exchanges 4,5,6. This is because in exchange 2,3 cultures were microwaved in 15 tubes for 3 minutes and in exchange 4,5,6 cultures were microwaved in 4 tubes for only 1 minute.
Tidy
#Columns: exchange & evaporation
ds_for_evaporation = gather(ds_for_evaporation,
key = exchange,
value = evaporation,
water_add_after_t2:water_add_after_t6)
ds_for_evaporation[ds_for_evaporation == "water_add_after_t2"] = "2"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t3"] = "3"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t4"] = "4"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t5"] = "5"
ds_for_evaporation[ds_for_evaporation == "water_add_after_t6"] = "6"
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] / 2 #This is because exchange contained the topping up of two exchanges
ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] = ds_for_evaporation$evaporation[ds_for_evaporation$exchange == 2] + 2 #We need to add 2 ml to the evaporation that happened at the exchange events 1 and 2. This is because we already added 1 ml of water at exchange 1 and 1 ml of water at exchange 2.
#Column: nr_of_tubes_in_rack
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 1] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 2] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 3] = 15
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 4] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 5] = 4
ds_for_evaporation$nr_of_tubes_in_rack[ds_for_evaporation$exchange == 6] = 4Plot
ds_for_evaporation %>%
filter(disturbance == disturbance) %>%
ggplot(aes(x = as.character(nr_of_tubes_in_rack),
y = evaporation)) +
geom_boxplot() +
labs(x = "Number of tubes in rack",
y = "Evaporation (ml)")
ds_for_evaporation %>%
filter(disturbance == disturbance) %>%
ggplot(aes(x = as.character(patch_size),
y = evaporation)) +
geom_boxplot() +
labs(x = "Patch size",
y = "Evaporation (ml)")
ds_for_evaporation %>%
filter(disturbance == disturbance) %>%
ggplot(aes(x = as.character(day),
y = evaporation)) +
geom_boxplot() +
labs(x = "Day",
y = "Evaporation (ml)")
ds_for_evaporation %>%
filter(disturbance == disturbance) %>%
ggplot(aes(x = disturbance,
y = evaporation)) +
geom_boxplot() +
labs(x = "Disturbance",
y = "Evaporation (ml)")
It seems like there is no real difference across time, disturbance, or patch type. However, we could also run a mixed effect model to show that they do not.
Mixed effect model
It gives me the following error:
- Error in fn(nM$xeval()) : Downdated VtV is not positive definite
mixed.model = lmer(evaporation ~
patch_size * disturbance * exchange +
(exchange | culture_ID),
data = ds_for_evaporation,
REML = FALSE,
control = lmerControl (optimizer = "Nelder_Mead"))
null.model = lm(evaporation ~
1,
data = ds_for_evaporation)
anova(mixed.model, null.model)Body size
Aim
Data
Experimental cultures
culture_info = read.csv(here("data", "PatchSizePilot_culture_info.csv"), header = TRUE)
load(here("data", "morphology", "t0.RData"));t0 = morph_mvt
load(here("data", "morphology", "t1.RData"));t1 = morph_mvt
load(here("data", "morphology", "t2.RData"));t2 = morph_mvt
load(here("data", "morphology", "t3.RData"));t3 = morph_mvt
load(here("data", "morphology", "t4.RData"));t4 = morph_mvt
load(here("data", "morphology", "t5.RData"));t5 = morph_mvt
load(here("data", "morphology", "t6.RData"));t6 = morph_mvt
load(here("data", "morphology", "t7.RData"));t7 = morph_mvt
rm(morph_mvt)
datatable(culture_info[,1:10],
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))Body size data-set
### --- Tidy t0 - t7 data-sets --- ###
#Column: time
t0$time = NA
t1$time = NA
#Column: replicate_video
t0$replicate_video[t0$file == "sample_00001"] = 1
t0$replicate_video[t0$file == "sample_00002"] = 2
t0$replicate_video[t0$file == "sample_00003"] = 3
t0$replicate_video[t0$file == "sample_00004"] = 4
t0$replicate_video[t0$file == "sample_00005"] = 5
t0$replicate_video[t0$file == "sample_00006"] = 6
t0$replicate_video[t0$file == "sample_00007"] = 7
t0$replicate_video[t0$file == "sample_00008"] = 8
t0$replicate_video[t0$file == "sample_00009"] = 9
t0$replicate_video[t0$file == "sample_00010"] = 10
t0$replicate_video[t0$file == "sample_00011"] = 11
t0$replicate_video[t0$file == "sample_00012"] = 12
t1$replicate_video = 1 #In t1 I took only 1 video/culture
t2$replicate_video = 1 #In t2 I took only 1 video/culture
t3$replicate_video = 1 #In t3 I took only 1 video/culture
t4$replicate_video = 1 #In t4 I took only 1 video/culture
t5$replicate_video = 1 #In t5 I took only 1 video/culture
t6 = t6 %>% rename(replicate_video = replicate)
t7 = t7 %>% rename(replicate_video = replicate)
### --- Create ds_body_size dataset --- ###
long_t0 = t0 %>% slice(rep(1:n(), max(culture_info$culture_ID)))
ID_vector = NULL
ID_vector_elongating = NULL
for (ID in 1:max(culture_info$culture_ID)){
ID_vector = rep(ID, times = nrow(t0))
ID_vector_elongating = c(ID_vector_elongating, ID_vector)
}
long_t0$culture_ID = ID_vector_elongating
t0 = merge(culture_info,long_t0, by="culture_ID"); rm(long_t0)
t1 = merge(culture_info,t1,by="culture_ID")
t2 = merge(culture_info,t2,by="culture_ID")
t3 = merge(culture_info,t3,by="culture_ID")
t4 = merge(culture_info,t4,by="culture_ID")
t5 = merge(culture_info,t5,by="culture_ID")
t6 = merge(culture_info,t6,by="culture_ID")
t7 = merge(culture_info,t7,by="culture_ID")
ds_body_size = rbind(t0, t1, t2, t3, t4, t5, t6, t7); rm(t0, t1, t2, t3, t4, t5, t6, t7)
### --- Tidy ds_body_size data-set --- ###
#Column: day
ds_body_size$day = ds_body_size$time_point;
ds_body_size$day[ds_body_size$day=="t0"] = "0"
ds_body_size$day[ds_body_size$day=="t1"] = "4"
ds_body_size$day[ds_body_size$day=="t2"] = "8"
ds_body_size$day[ds_body_size$day=="t3"] = "12"
ds_body_size$day[ds_body_size$day=="t4"] = "16"
ds_body_size$day[ds_body_size$day=="t5"] = "20"
ds_body_size$day[ds_body_size$day=="t6"] = "24"
ds_body_size$day[ds_body_size$day=="t7"] = "28"
ds_body_size$day = as.numeric(ds_body_size$day)
#Column: time point
ds_body_size$time_point[ds_body_size$time_point=="t0"] = 0
ds_body_size$time_point[ds_body_size$time_point=="t1"] = 1
ds_body_size$time_point[ds_body_size$time_point=="t2"] = 2
ds_body_size$time_point[ds_body_size$time_point=="t3"] = 3
ds_body_size$time_point[ds_body_size$time_point=="t4"] = 4
ds_body_size$time_point[ds_body_size$time_point=="t5"] = 5
ds_body_size$time_point[ds_body_size$time_point=="t6"] = 6
ds_body_size$time_point[ds_body_size$time_point=="t7"] = 7
ds_body_size$time_point = as.character(ds_body_size$time_point)
#Column: eco_metaeco_type
ds_body_size$eco_metaeco_type = factor(ds_body_size$eco_metaeco_type,
levels=c('S', 'S (S_S)', 'S (S_L)', 'M', 'M (M_M)', 'L', 'L (L_L)', 'L (S_L)'))
#Select useful columns
ds_body_size = ds_body_size %>%
select(culture_ID,
patch_size,
disturbance,
metaecosystem_type,
mean_area,
replicate_video,
day,
metaecosystem,
system_nr,
eco_metaeco_type)
#Reorder columns
ds_body_size = ds_body_size[, c("culture_ID",
"system_nr",
"disturbance",
"day",
"patch_size",
"metaecosystem",
"metaecosystem_type",
"eco_metaeco_type",
"replicate_video",
"mean_area")]
datatable(ds_body_size,
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))## Warning in instance$preRenderHook(instance): It seems your data is too big
## for client-side DataTables. You may consider server-side processing: https://
## rstudio.github.io/DT/server.htmlSize classes data-set
I am here creating 12 size classes as in Jacquet, Gounand, and Altermatt (2020). However, for some reason it seems like our body size classes are really different.
#### --- PARAMETERS & INITIALISATION --- ###
nr_of_size_classes = 12
largest_size = max(ds_body_size$mean_area)
size_class_width = largest_size/nr_of_size_classes
size_class = NULL
### --- CREATE DATASET --- ###
size_class_boundaries = seq(0, largest_size, by = size_class_width)
for (class in 1:nr_of_size_classes){
bin_lower_limit = size_class_boundaries[class]
bin_upper_limit = size_class_boundaries[class+1]
size_input = (size_class_boundaries[class] + size_class_boundaries[class + 1])/2
size_class[[class]] = ds_body_size%>%
filter(bin_lower_limit <= mean_area) %>%
filter(mean_area <= bin_upper_limit) %>%
group_by(culture_ID,
system_nr,
disturbance,
day,
patch_size,
metaecosystem,
metaecosystem_type,
eco_metaeco_type,
replicate_video) %>% #Group by video
summarise(mean_abundance_across_videos = n()) %>%
group_by(culture_ID,
system_nr,
disturbance,
day,
patch_size,
metaecosystem,
metaecosystem_type,
eco_metaeco_type) %>% #Group by ID
summarise(abundance = mean(mean_abundance_across_videos)) %>%
mutate(log_abundance = log(abundance)) %>%
mutate(size_class = class) %>%
mutate(size = size_input) %>%
mutate(log_size = log(size))
}
ds_classes = rbind(size_class[[1]], size_class[[2]], size_class[[3]], size_class[[4]],
size_class[[5]], size_class[[6]], size_class[[7]], size_class[[8]],
size_class[[9]], size_class[[10]], size_class[[11]], size_class[[12]],)
datatable(ds_classes,
rownames = FALSE,
options = list(scrollX = TRUE),
filter = list(position = 'top',
clear = FALSE))Plot
Comparison of different ecosystem types across time
#Trying out gganimate, but I can't seem to manage to install transformr packaget
p = list()
n = 0
first_level = c("isolated small", "isolated small", "isolated large", "isolated large")
second_level = c("small connected to small", "small connected to small", "large connected to large", "large connected to large")
third_level = c("small connected to large", "small connected to large", "large connected to small", "large connected to small")
for (patch_size_input in c("S", "L")){
for(disturbance_input in c("low", "high")){
n = n + 1
title = paste0(patch_size_input,
' patches, Disturbance = ',
disturbance_input,
', Day: {round(frame_time, digits = 0)}')
p[[n]] <- ds_classes %>%
filter(disturbance == disturbance_input) %>%
filter(patch_size == patch_size_input) %>%
ggplot(aes(x = log_size,
y = log_abundance,
group = interaction(log_size, eco_metaeco_type),
color = eco_metaeco_type)) +
geom_point(stat = "summary", fun = "mean") +
geom_line(stat = "summary", fun = "mean", aes(group=eco_metaeco_type)) +
scale_color_discrete(labels = c(first_level[n],
second_level[n],
third_level[n])) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
legend.position = c(.95, .95),
legend.justification = c("right", "top"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)) +
labs(title = title,
x = 'Log size (μm2)',
y = 'Log abundance + 1 (indiv/μm2)',
color = "") +
transition_time(day) +
ease_aes('linear')
animate(p[[n]],
duration = 10,
fps = 25,
width = 500,
height = 500,
renderer = gifski_renderer())
anim_save(here("gifs",
paste0("transition_day_",
patch_size_input,"_",
disturbance_input,
".gif")))
}
}
- Next step: making plots slower (the package stopped working because I broke some dependencies)
Figures
Regional biomass under low disturbance in meta-ecosystems of the same total area, but exhibiting two local patches of the same size (red, medium-medium meta-ecosystem) or different size (blue, small-large meta-ecosystem). Points represent the mean, error bars represent the standard deviation. See the Appendix for equivalent figure of the high disturbance treatment.
Local biomass under low disturbance in patches connected to a patch of the same size (red, small patches of small-small meta-ecosystems) or to a patch of larger size (blue, small patches of small-large meta-ecosystems). The local biomass has been transformed as the natural logarithm of the reponse ratio as follows: ln(RR) = ln((bioarea + 1) / mean bioarea small isolated patches under low disturbance).
Local biomass under low disturbance in patches connected to a patch of the same size (red, large patches of large-large meta-ecosystems) or to a patch of larger size (blue, large patches of small-large meta-ecosystems). The local biomass has been transformed as the natural logarithm of the reponse ratio as follows: ln(RR) = ln((bioarea + 1) / mean bioarea large isolated patches under low disturbance).
Tests
Evaporation when microwaving 15 falcon tubes at the time
evaporation.test = read.csv(here("data", "evaporation_test","evaporation_test_right.csv"), header = TRUE)
evaporation.test %>%
ggplot(aes (x = as.character(water_pipetted),
y = weight_water_evaporated,
group = interaction(water_pipetted, as.character(rack)),
fill = as.character(rack))) +
geom_boxplot() +
labs(x = "Water volume (ml)" ,
y = "Evaporation (g)",
fill = "Rack replicate")
Evaporation when microwaving 5 tubes with 10 filled or empty tubes
evaporation.test = read.csv(here("data", "evaporation_test", "evaporation_test_fill_nofill.csv"), header = TRUE)
evaporation.test %>%
ggplot(aes (x = all_tubes_water,
y = weight_water_evaporated)) +
geom_boxplot() +
labs(x = "Water in the other 10 tubes" ,
y = "Evaporation (g)",
caption = "When all tubes were filled, they were filled with 6.75 ml of deionised water.")
Running time
## Time difference of 37.18174 secsOther
Mixed effects models
To build the mixed effect models we will use the R package lme4. See page 6 of this PDF to know more about the syntaxis of this package and this link for the interaction syntaxis.
To do model diagnostics of mixed effect models, I’m going to look at the following two plots (as suggested by Zuur et al. (2009), page 487):
Quantile-quantile plots (
plot(mixed_model))Partial residual plots (
qqnorm(resid(mixed_model)))
The effect size of the explaining variables is calculated in the mixed effect models as marginal and conditional r squared. The marginal r squared is how much variance is explained by the fixed effects. The conditional r squared is how much variance is explained by the fixed and the random effects. The marginal and conditional r squared are calculated using the package
MuMIn. The computation is based on the methods of Nakagawa, Johnson, and Schielzeth (2017). For the coding and interpretation of these r squared check the documentation for the r.squaredGLMM functionTime can be included as a fixed or random effect. Time can be included as a random effect if the different data points are non independent from each other (e.g., seasons). However, because the biomass in our experiment was following a temporal trend, the different time points show autocorrelation. In other words, t2 is more similar to t3 than t4 and so on. This is why we decided to include time as a fixed effect. For an excellent discussion on this topic see this blog post.
Modeling choices
I am going to select the best model according to AIC. Halsey (2019) suggests this approach instead of p values. P-values are not a reliable way of choosing a model because:
My sample size is small, producing larger p-values
P-values are really variable, creating many false positives and negatives (e.g., if p=0.05 there is a 1 in 3 chance that it’s a false positive)
To study the local biomass how it changes across treatments, we could have made three different models between the three combinations of small patches. However, that might be confusing to interpret the results. We decided instead to use an effect size where we control is the isolated small patch. At the beginning we thought to use the natural logarithm of the response ratio (lnRR). The problem, however, is that some bioarea values were 0. We were thinking to add 1 to all null values, but according to Rosenberg, Rothstein, and Gurevitch (2013), such practice inflates effect sizes. Because of this, I looked into other types of effect size. I found that the most common and preferred metric in use today is known as Hedge’s d (a.k.a. Hedge’s g) (Hedges, Larry V. and Olkin (1985) ). It is calculated as the difference in mean between treatment and control divided by the standard deviation of the pooled data. Another measure would be Cohen’s d, but it underperforms with sample sizes that are lower than 20 (StatisticsHowTo). I can easily calculate the Hedge’s d using the r package
effsize.
Same thing for the large patches.
Bibliography
Halsey, Lewis G. 2019. “The reign of the
p-value is over: What alternative analyses could we employ to fill the
power vacuum?” Biology Letters 15 (5). https://doi.org/10.1098/rsbl.2019.0174.
Hedges, Larry V., and Ingram Olkin. 1985. Statistical Methods for Meta-Analysis.
Jacquet, Claire, Isabelle Gounand, and Florian Altermatt. 2020.
“How pulse disturbances shape size-abundance
pyramids.” Ecology Letters 23 (6): 1014–23. https://doi.org/10.1111/ele.13508.
Nakagawa, Shinichi, Paul C. D. Johnson, and Holger Schielzeth. 2017.
“The coefficient of determination R2 and
intra-class correlation coefficient from generalized linear
mixed-effects models revisited and expanded.” Journal
of the Royal Society Interface 14 (134). https://doi.org/10.1098/rsif.2017.0213.
Rosenberg, Michael S., Hannah R. Rothstein, and Jessica Gurevitch. 2013.
“Effect sizes: Conventional choices and
calculations.” Handbook of Meta-Analysis in Ecology
and Evolution, 61–71. https://doi.org/10.23943/princeton/9780691137285.003.0006.
Zuur, Alain F., Elena N. Ieno, Neil Walker, Anatoly A. Saveliev, and
Graham M. Smith. 2009. Mixed effects models and
extensions in ecology with R. Vol. 36. Statistics for
Biology and Health. New York, NY: Springer New York. https://doi.org/10.1007/978-0-387-87458-6.